はじめに

実習の目的

生物個体群における種内・種間の形態変異パターンを解析し、その理解と方法の習得を目指す。各自に与えられた二枚貝類のサンプルについて、観察、計測、集計、解析などを行い、その結果をレポートにまとめる。

使用ツール

本実習では、表計算ソフトの Microsoft Excel と統計解析用フリーソフト R を使用する。 Excel は、データの編集や作図を直感的かつ視覚的に行える一般的な表計算ソフトである。多数の関数が利用可能で、小規模なデータの集計や基本的な作図に優れている。一方で、複雑なグラフの作成や高度な統計解析には向いていない。 R は、先端的な統計解析まで対応可能な汎用性の高いソフトウェアである。行列操作やコマンドラインでの操作に慣れが必要だが、データの集計や作図も含めて幅広い解析に対応できる。本実習ではExcelとRを適宜使い分けながら、基礎的なデータ解析手法を習得する。

生物測定学・生物統計学の意義

生物測定学（biometry） および 生物統計学（biostatistics） は、生物のさまざまな形質や特性を定量的に測定・解析する学問である。生物学の1分野というよりも、生物科学全般において定量データから客観的な結論を導くための普遍的な手段と位置づけられる。この実習では、その初歩的な考え方と方法を体験的に学ぶ。これらの手法は、科学的なデータ処理や実験計画、さらには製品比較や品質管理など、実社会での応用にも役立つ。

個体群とサンプル解析の意義

個体群（population）は、生物が生活するうえでの基本単位である。遺伝、生態、分類、進化を考察する際にも、同じ場所で繁殖社会を形成している自然の個体群の性質を正しく理解することが出発点となる。自然界の個体群すべてを調査することは現実的に困難であるため、個体群からランダムに抽出されたサンプル（sample）を解析し、個体群全体の性質を統計的に推定する。

本実習で扱うサンプル

本来、サンプルは各自が野外で生態を観察しながら、母集団からランダムに採集すべきものである。しかし、時間的制約のため、あらかじめ準備された旧金沢研究室所有の二枚貝類標本を使用する。これらの標本は、同一の産地からランダムに採集されたと考えられる。他人のサンプルと混同したり、持ち帰ったりしてはならない。各自が異なるサンプルを担当し、それぞれ独自に計測データを得て解析を行う。課題の結果は、サンプルごとに多少異なる場合がある。基本的には、本テキストの指示に従い、種内・種間の形態変異パターンを分析する。


実習実施予定

予定

  • 第1部 サンプルの前処理、形質の定量化（計測、集計、計測データ入力）
  • 第2部 データの集計
  • 第3部 種内変異の解析
  • 第4部 種間変異の解析

用いる道具と事前準備

  • 電子ノギス（一部の大型種については別途大型のノギスを貸し出す）
  • 集計用紙
  • タックラベル（貝殻の左右、個体識別用）
  • 補助板（貝殻の厚みの計測用）
  • コンピュータ（データ解析用、Windows PCの使用を推奨）

事前に最新バージョンのExcelが動作する様にしておく事（JINDAIアカウントに紐づけられるOffice365を利用して無料でインストールする事ができる）。また、Rを事前にダウンロード、インストールして、使えるようにしておくこと。Rのダウンロードは以下のリンクより行うこと。Rのインストールは、全てデフォルトの設定で行う。

R（Windows版）のインストール手順

基本的にはデフォルト設定で、全て「次へ」でOK。 R_install



第1部 サンプルの前処理、形質の定量化

形態変異は何らかの形質（character）を対象として解析する。形質にはさまざまなものがあるが、年齢の不揃いなサンプルの変異の解析には成長にともなって変化しにくいものを選ぶ必要がある。例えば、殻の大きさや重量そのものは、成長にともなって次第に増加するから、変異の解析として集計しても意味がないことが多い。これに対して、長さと高さの単純比（simple ratio）や相同的な器官の数は、成長にともなって変化しにくいと考えられるため、個体変異を表す形質として扱うことができる。二枚貝の外表にしばしば見られる放射肋(ほうしゃろく)の数は、もし分岐や挿入がなければ、意味のある形質となる。これらは1変量形質（univariate character）または2変量の単純比を1変量形質として扱うものである。この他に2変量形質（bivariate character）、多変量形質（multivariate character）を解析することもある。

1-1. サンプルの前処理

1-1-1. 種名の確認

大部分のサンプルには種名が一応与えてあるが、異なった種が混合していないか各自で確かめること。複数種が混合している可能性がある場合には、種名が正しいか計測終了までに随時インターネットでの検索やスタッフへの確認などで確かめておくこと。また、種名カードが箱に入っていない場合には、インターネット等で種名を検索し、種名の特定を行うこと。属名や種名は図鑑の著者の見解によって異なることがあり、変更されている事もある。

1-1-2. 生活様式によるグループ分け

二枚貝はその生活様式を大まかに次の3種類に分類することが出来る：①堆積物に潜るタイプ、②堆積物の表面で生活するタイプ、③基盤に付着するタイプ。ここでは、①の生活様式を持つ種をAグループ、②と③の生活様式を持つ種をBグループとして、自分のサンプルがAとBのどちらの生活様式グループに属するかを以下の図Aを参考にして判定する。

AB_grouping 図A 生活様式AグループとBグループの判別方法。

1-1-3. 二枚貝の方位の決定

二枚貝には、前後対称の種はない。ハマグリのようにほぼ左右対称の種が多いが、ホタテガイのように左右非対称の種も少なくない。以下の識別点を参考にして、二枚貝の方位（前後左右）を認識し、左右の殻を区別する。分からない場合にはスタッフに確認する。

  1. 通常の ハマグリ型の種では、殻頂は中央よりも前方に位置し（例外あり）、後縁が裁断状になる。
  2. 套線湾入(とうせんわんにゅう)（水管を格納するスペース）があれば必ず後側にできる。
  3. 一般に靭帯(じんたい)は殻頂の後側にできる（例外あり）。
  4. ホタテガイ類では、前方の耳状部がより突出し、前側の閉殻筋(へいかくきん)が退化する。また、一般に左殻の方が色が濃い。
  5. 殻頂を上に、殻の前部を遠方に向けた時、右側の殻が右殻、左側の殻が左殻である。
  6. 前後、左右ともほぼ対称の種もあるが、注意すれば見分けられる。

1-1-4. サンプルの取り扱い方

  1. 扱うサンプルの個体数（sample size）は多いほどよい（誤差が小さくなる）。ここでは20個体（貝殻の総数がそれ以下の場合は全個体）を原則として、時間がある場合には25個体を計測する。種によっては個体数が多いので、その中から20個体をランダムに抽出して用いる
  2. 各個体の外面または内面にタックラベルを貼って番号を書く。貝殻に凹凸が多い場合、外側だと剝がれやすいので、注意すること。貝殻の左右と標本番号を併せて、右殻はR1、R2、R3…、左殻はL1、L2、L3…、のように記す。
  3. 貝の形状に合わせて、測定の基準となる「基線」を定める。前後の閉殻筋痕の中心を結ぶ線や、殻頂部分の直線などが基線の候補となり得る。

1-2. 計測

長さ（L）、高さ（H）、厚さ（T）、右殻厚さ（TR）、左殻厚さ（TL）をノギスで計測する。作業を開始する前に、貝の基線とノギスの操作方法を確認しておくこと。ノギスの操作・計測に1つでも誤り（表示値の読み違い、転記ミスなど）があると、以後の作業がすべて無駄になる可能性があるので十分注意すること。先に定めた基線に平行または直交する方向にLとHを測定する（図B）。合弁の個体は、まず殻を合わせた状態でTを計測し、次に右殻と左殻を別々に計測してTR、TLとするのがよい。T=TR+TLにならないことが多いが、それは差し支えない。二枚貝類の殻の厚さはノギスでは測りにくいが、補助板を当てるなど各自で工夫すること。補助板を用いる場合には、補助板の厚みを測定し、後で補助版込みの測定値から差し引くことで補正しておくこと。LとHは左右共通の値として、左右いずれかの殻に決めて測定を行う。計測は独立に2回行って数値を個体ごとに集計用紙に記録する。5個体ずつを一つのブロックにし、5個体についての1回目の計測が終わったら、同じ5個体について2回目の計測を行い、2回の計測が終了した後に、次の5個体のブロックへと進めるようにする。使用しているノギスの精度に応じて、小数点以下１桁もしくは２桁で記録を行う。

各自に割り振られた貝について、基線の位置、L、H、T、TR、TLの測定位置を写真もしくはイラストにして提出すること（課題１, Wordファイル）。

measurement1 図B 計測部位の例。L：殻長、H：殻高、TR：右殻厚


measurement2 図C 各部部位の計測方法


集計用紙の記入例： measurement3



第2部 データの集計

2-1. 測定値の入力と平均値の算出

  1. シート「元データ」には個体番号1-20（20個体分に満たない場合には、ある分だけ）のそれぞれの個体について、各パラメータの測定値２回分を入力する。入力の単位はmmとする。
  2. 2回の反復測定の平均値を算出し、「解析用データ」シートに入力する。

測定には誤差がつきものであり、多数回の測定を行うことで誤差の影響を低減することができるが大きな労力が必要となる。この実習では2回の測定値を平均することで簡易的に誤差を低減させた値をもとに形態解析を行なう。2回の測定値に大きな違いがある場合には、貝の取違いや、測定部位に誤りがある可能性が非常に高いので、確認をすること。測定に用いるノギスによって、測定精度は1/10mm～１/100mm程度である。使用したノギスの制度に合わせて、Excelのround関数とaverage関数を組み合わせて平均値を1/10mmもしくは1/100mm精度として以降の処理に用いる。

平均値の算出（Excelを使用) 動画1

「元データ」シートからデータをコピー＆ペーストして、平均値を格納したい表の隣に配置する。その後、average関数を用いて、対応するセル間の平均値を計算する。また、平均値を四捨五入して、小数点以下2桁で表現するために、round関数を使用する。平均化と四捨五入の処理を同時に行うためには、「round(average(○○:○○),2)」（ただし、○○部分は対応するセル）と入力すればよい。関数の使用時には、最初に「＝」の入力を忘れないこと。

average　



第3部 種内変異の解析

個体群中の個体は、一般的にそれぞれ異なった遺伝的性質を持ち、少しずつ異なった環境の影響を受けている。そのために、形態的に完全に同一な個体は存在しない。たとえ成長段階が揃っていても、形態は多くの遺伝的要因と環境的要因の相互作用によって支配され、変異（variation）が生じている。単一または少数の遺伝的要因が形態に特別に強い影響を与えるときには不連続変異（例えば、性染色体による雌雄の違いや対立遺伝子による多型現象）も生ずる。単一の個体群中にみられる変異を個体変異（individual variation）といい、同一種の異なる個体群間で見られる変異を地理的変異（geographic variation）という。この実習で取り扱う貝殻は、基本的に単一の個体群から採集されていると仮定し、形態の種内変異は個体変異であると考える事とする。通常の連続的な個体変異には、遺伝的変異（突然変異に起因する）と環境的変異（遺伝しない変異）の両方が含まれているはずであるが、両者を量的に識別することはこの実習では困難であるため、ここでは両者を区別しない。

3-1. 変量の比と代表値の算出

HとL、T、TR、TLの値を2つずつ組み合わせた比（L/H、T/H、TR/H、TL/H） を求める。さらに各形質（2変量の比）について、次の代表値を計算せよ（課題2、Excelファイル）。

  1. 標本数（\(N\)）:\(N\)
  2. 平均値（\(\bar{x}\)): \(\sum \frac{x_i}{N}\)
  3. 標準偏差（\(s\)）:\(\sqrt{ \frac{1}{N - 1} \sum (x_i - \bar{x})^2 }\)
  4. 平均値の標準誤差（\(s_{\bar{x}}\) ）:\(\frac{s}{\sqrt{N}}\)
  5. 変動係数（\(V\)）:\(\frac{100s}{\bar{x}}\)
  6. 最小値（\(X_{min}\)）
  7. 最大値（\(X_{max}\)）

\(N\)、\(\bar{x}\)、\(s\)はサンプルからの母集団（個体群）の性質の推定や、他のサンプルとの比較 を行うときに最も基本となる代表値である。

各代表値の計算方法（Excelを使用） 動画2

  • 標本数（\(N\)）:データ数をカウントする。今回は不要だが、count関数を使うこともできる。
  • 平均値（\(\bar{x}\)):：average関数を用いて計算する。
  • 標準偏差（\(s\)）：stdev.s関数を使用して計算する。
  • 平均値の標準誤差（\(s_{\bar{x}}\) ）：sqrt関数を使い、s/sqrt(N)として計算する。
  • 変動係数（\(V\)）:平均値と標準偏差を使って計算する。
  • 最小値（\(X_{min}\)）：min関数を使って計算する。
  • 最大値（\(X_{max}\)）：max関数を使って計算する。

Excel上での計算ミスは、ほとんどの場合入力ミスあるいは関数や数式の誤った使用によって生じる。また、表示桁数が多い事は計算上有利であるが、視認性が悪くなり異常なデータを認識しづらくなる。ここでは視認性を高めるため、小数点以下4桁程度に表示される様にしておく（右クリック→「セルの書式設定」→「数値」→「小数点以下の桁数」、ROUND関数とは異なり、データの実体は変更されない）。

ノギス測定時には明確な有効数字が存在するので、その後の計算値を示す際の有効数字はそれを基準に定めることになる。手計算や計算尺、電卓の時代とは異なり、計算そのものの煩雑さが問題とならない現代では、後の作業過程に充分以上の桁数を確保しておき、表示値を調整する場合が多い（厳密な有効数字とは対立するが、扱いやすい方法である）。今回は有効桁数を考慮せずに解析を行う。

3-2. データ分布の視覚化

生物測定学的研究では、データが正規分布（normal distribution）を示すこと（=サンプルが正規分布する母集団からランダムにサンプリングされていること）を前提に議論を進めることが多い。データが正規分布に近いということは、それが単一の種の個体群を代表している可能性が高いことを意味する。一般に、量的形質の分布は正規分布になりやすい。これは、量的形質が効果の小さい複数の遺伝子座（量的形質遺伝子座と呼ぶ）によって支配されていること、環境の影響を受けるとともに、遺伝子型と環境の相互作用が影響してるためと言われている。その一方で、質的形質を支配する遺伝子座は1つまたは少数で、環境の影響をほとんど受けない点で、量的形質と明確に異なっている。

正規分布は、多数の要因によって支配される生物の量的形質によく現れる分布型である。その認定は以後の統計解析を進める上で重要な前提になる。まずは、ヒストグラム（度数分布図、図D）を作って、サンプルの各形質が正規分布すると見なせるか視覚的に判定する。今回のように単純比のような連続変量を扱うときにはヒストグラムの階級の数を\(\sqrt{N}\)程度とする（Nはサンプル数）。階級を細かく分け過ぎると偶然性のためにデータの分布を認識しにくくなる。

図D ヒストグラムの例。

L/H、T/H、TR/H、TL/Hについてヒストグラムを作成せよ（課題３、Excelファイル）。

ヒストグラムの作成方法（ExcelとRを使用） 動画3

Rを起動して、以下のコードをコピー＆ペーストする。

data <- read.table("clipboard", header = TRUE)

＊Macの場合にはこちら：　data<-read.table(pipe(“pbpaste”),header=T)

次に、Excel上でL/H、T/H、TR/H、TL/Hの行列全体をコピーする（ただし、ペーストはしない）。

その後で、R上で先ほどのコードを実行（enterキーを押す）することで、Ｒ上でdataにデータ行列が格納される。データの格納に成功した場合には、以下のコードでデータ行列の最初の６行が表示される。

head(data)
##      L.H    T.H   TR.H   TL.H
## 1 1.0163 0.2562 0.1613 0.1300
## 2 0.9796 0.2447 0.1475 0.1101
## 3 0.9791 0.2334 0.1313 0.1089
## 4 0.9702 0.2432 0.1303 0.1126
## 5 1.0270 0.2617 0.1485 0.1376
## 6 0.9823 0.2619 0.1428 0.1288

R上では、列名の「/」は「.」に自動的に変換される。ヒストグラムはhist関数を用いて描画する。L/H、T/H、TR/H、TL/Hのヒストグラムの描画には、以下のコードを用いる。

# L/H(結果を例示)
hist(data$L.H, xlab="L/H", ylab="頻度", main="") 

# T/H
hist(data$T.H, xlab="T/H", ylab="頻度", main="") 
# TR/H
hist(data$TR.H, xlab="TR/H", ylab="頻度", main="") 
# TL/H 
hist(data$TL.H, xlab="TL/H", ylab="頻度", main="")

＊Macで漢字が文字化けする場合には、”頻度”を”frequency”に置き換える。

コードを実行して得られたヒストグラムのグラフ上で右クリックをして、グラフを「メタファイルにコピー」する。＊Macの場合には、やや方法が異なる。コピーの方法が分からない場合には、スクリーンショットを貼り付けても良い。

その後、Excel上の「ヒストグラム」シート上の指定の位置に貼り付ける。

ヒストグラムを眺めて、正現分布になりそうかどうか推測せよ（課題4、Wordファイル）。また、サンプルの中に他と大きく異なった個体がないか、もしあればその個体を他と同一の種と見なせるか、改めてその個体に戻って他の形質と関連があるか検討せよ。ここで計測や計算の誤りが見つかるかも知れない。

hist関数に引数nclassを追加して、階級数を変化させたときのグラフも2つ追加し、階級数がヒストグラムの見た目に与える影響についても考察せよ（追加課題1と2、Excelファイル）。この時、ヒストグラム描画に使用するデータはL/H、T/H、TR/H、TL/Hのいずれでも構わない。また、指定する階級数も自由に設定して構わない。ただし、使用したデータと階級数をExcelシートに記載すること。 ＊生物科学実験IIの履修者は必須の課題とする。



3-3. 正規分布とのずれの視覚化

正規分布はヒストグラムを見て直観的に判断することもあるが、もう少し理論的に厳密に確認する必要がある。特に個体数が比較的少ない場合には正規分布になるか判定しにくいことが多い。一つの方法として考えられるのは、母集団が正規分布であると仮定したときの理論的な度数分布と実際の観測度数分布を総合的に比較することである。ヒストグラムに観測値から算出した平均値と標準偏差を持つ正規分布を上書きし、正規分布とのずれを視覚化する（課題5、Excelファイル）。

ヒストグラムへの正規分布の上書き（Rを使用） 動画4

先ほどと同様に、Ｒ上でdataにデータ行列を格納する。ヒストグラムはhist関数を用いて描画し、正規分布はcurve関数を用いる。描画には、以下のコードを用いる。

# L/H(結果を例示)
hist(data$L.H,xlab="L/H",ylab="確率密度",main="",freq=F) 
curve(dnorm(x,mean(data$L.H),sd(data$L.H)),add=TRUE,col="red")

# T/H
hist(data$T.H,xlab="T/H",ylab="確率密度",main="", freq=F) 
curve(dnorm(x,mean(data$T.H),sd(data$T.H)),add=TRUE,col="red")
# TR/H
hist(data$TR.H,xlab="TR/H",ylab="確率密度",main="", freq=F) 
curve(dnorm(x,mean(data$TR.H),sd(data$TR.H)),add=TRUE,col="red")
# TL/H 
hist(data$TL.H,xlab="TL/H",ylab="確率密度",main="", freq=F) 
curve(dnorm(x,mean(data$TL.H),sd(data$TL.H)),add=TRUE,col="red")

＊確率密度とは、確率変数の取りうる値に対して、その発生する確率を表したもの。確率変数の範囲全体での積分値が1となる。

＊Macで漢字が文字化けする場合には、”確率密度”を”density”に置き換える。

コードを実行して得られたヒストグラムのグラフ上で右クリックをして、グラフを「メタファイルにコピー」し、Excel上の「正規性」シート上の指定の位置に貼り付ける。＊Macの場合には、やや方法が異なる。



3-4. 検定とは？

一般に、統計学の検定は次のような考え方で行われる（どの検定でも共通）。

  1. 帰無仮説（\(H_0\)）を立てる
  • 「母集団は正規分布する」や「2つの母集団の平均値には差がない」のような唯一の状態（\(H_0\)）を仮定する。逆に、「正規分布でない」や「平均値に差がある」のような仮説には無限の可能性があり、帰無仮説には適さない。
  1. \(H_0\)が成り立つ確率（P値）を求める
  • 適切な統計手法に従って、データが帰無仮説のもとで得られる確率を計算する。
  1. P値に基づいて判断する
  • P < 0.05 の場合：\(H_0\)の正当性が疑われる ⇒ \(H_0\)を棄却し、対立仮説\(H_1\)を採択する
  • P ≥ 0.05 の場合：\(H_0\)の正当性を疑う根拠がない ⇒ \(H_0\)を採択する（ただし、H₀が正しいと断定するわけではない）

3-5. 分布の正規性の検定（KS検定）

旧ソビエト連邦の数学者 Andrey Nikolaevich KolmogorovとNikolai Vasilyevich Smirnov によって開発されたコルモゴロフ–スミルノフ検定（KS検定）は適合度検定の一つで、二つの頻度分布の差を検定することができる。正規性の検定に用いる場合には、一つは観測分布、もう一つは期待分布（正規分布を仮定した場合の期待値）を用いる。観測値の累積確率分布と正規分布の累積確率の差の絶対値の最大値Dを検定統計量とし検定を行う（図E）。

図E 観測値の累積確率分布と正規分布の累積確率の差の絶対値

KS検定によって正規性の検定を行い、P値が0.05を下回った場合には、帰無仮説（＝観測値は正規分布する）を棄却して、正規分布からずれていると結論する。L/H、T/H、TR/H、TL/Hについて、データの分布が正規分布からずれているかどうかを、KS検定によって検定し、P値をエクセルファイルに記入せよ（課題6、Excelファイル）。また、L/H、T/H、TR/H、TL/Hについて、KS検定の結果に基づいて、正規分布からのずれについての解釈と、その根拠を記述せよ（課題7、Wordファイル）。

KS検定を用いた正規性の検定（Rを使用） 動画5

先ほどと同様に、Ｒ上でdataにデータ行列を格納する。KS検定には、以下のコードを用いる。

# L/H(結果を例示)
ks.test(data$L.H, "pnorm", mean=mean(data$L.H), sd=sd(data$L.H))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  data$L.H
## D = 0.15955, p-value = 0.498
## alternative hypothesis: two-sided
#T/H
ks.test(data$T.H, "pnorm", mean=mean(data$T.H), sd=sd(data$T.H))
#TR/H
ks.test(data$TR.H, "pnorm", mean=mean(data$TR.H), sd=sd(data$TR.H))
#TL/H
ks.test(data$TL.H, "pnorm", mean=mean(data$TL.H), sd=sd(data$TL.H))

得られたp-valueの値をコピーして、「正規性」シートの該当箇所（「KS検定　P=」の部分に貼り付ける。

＊データ中に同じ数値（同順位の値）が存在していると、” ties should not be present for the one-sample Kolmogorov-Smirnov test“というエラーメッセージが表示されるが、ここでは気にせずに検定を行う。



3-6. 平均値の差の検定（t検定）

2つのサンプルのある量的形質の平均値には一般に多少の差がある。その差が意味のあるものかどうか（つまり2つの母集団での平均値が等しいかどうか）を統計学的に明らかにするために利用されるのがt検定である。ここでは、実際にt検定を行って左右の殻の膨らみ（TR/HとTL/H）の平均値に有意の差があるかどうか（左右の殻の膨らみに差があるかどうか）を統計学的に判定する。t検定は個体数が少なくても適用できるが、2つのサンプルの変量が共に正規分布することが前提となるので、事前にKS検定で正規分布を確かめておくこと。ただし、今回は実習であるから「正規分布ではない」という結果が得られた場合にもt検定を行うこと（その旨はレポートに明記しておくこと）。

2つのサンプルの平均値を\(\bar{x}_1\)、\(\bar{x}_2\)、標準偏差を\(s_1\)、\(s_2\)、個体数を\(N_1\)、\(N_2\)とすると、次の式で\(t\)値が求められる。

\[ t = \frac{(\bar{x}_1 - \bar{x}_2) \sqrt{\frac{N_1 N_2}{N_1 + N_2}}} {\sqrt{ \frac{(N_1 - 1) s_1^2 + (N_2 - 1) s_2^2}{N_1 + N_2 - 2} }} \] この\(t\)値から、帰無仮説（\(H_0\):\(μ_1\)=\(μ_2\)）が成り立つ確率（P）を知ることができる。P≥0.05であれば、「平均値の差異が有意であるとはいえない」、P<0.05であれば「平均値の差異は有意である」という結論に達する。

TR/HとTL/Hについて、平均値±標準誤差の棒グラフを作成し、その差をt検定で検定せよ（課題8、Excelファイル）。また、その解釈を記述せよ（課題9、Wordファイル）。

棒グラフ（エラーバーつき）の作図 動画6

「変量の比と代表値」シートのTR/HとTL/Hの平均値を選択した状態で、Excelで「挿入」タブから、棒グラフを選択する。

出現したグラフエリアを右クリックし、「データの選択」を選び、「横（項目）軸ラベル（C）」の編集タブをクリック。



「軸ラベルの範囲（A）：」の上矢印をクリックし、TR/HとTL/Hを選択する。



横軸のラベルが設定できたら、誤差範囲の設定を行う。グラフエリアの右上の「＋」アイコン、その後「誤差範囲」をクリックする。



ここで表示されるエラーバーは誤差として不適切なので、どちらかのエラーバーを右クリックして、「誤差範囲の書式設定（E）」を選択する。



ここで表示される誤差範囲の中から、「ユーザー設定（C）」を選び、値を指定する。



正の誤差の値、負の誤差の値のいずれにもTR/HとTL/Hの標準誤差を指定する。



自分で算出した標準誤差を、以下のように選択する。



得られたグラフを「平均値の検定」シートの枠内にカット＆ペーストする。時間があれば、縦軸の追加や、枠線の削除など、グラフの見た目をきれいにする処置を行ってもよい。



t検定 動画7

これまでと同様に、Ｒ上でdataにデータ行列を格納する。t検定には、以下のコードを用いる。

t.test(data$TR.H, data$TL.H)
## 
##  Welch Two Sample t-test
## 
## data:  data$TR.H and data$TL.H
## t = 8.6032, df = 47.964, p-value = 2.759e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.03361548 0.05412052
## sample estimates:
## mean of x mean of y 
##  0.157680  0.113812

得られた結果から、p-valueの値をコピーして、「平均値の検定」シートの該当箇所（「t検定　P=」の部分）に貼り付ける。例えば、P=2.796e-11（e-11は、\(10^{-11}\)の意味）のような表記になっている場合でも、そのまま貼り付けて良い。出力で「p-value <」と表示されている場合には、エクセルファイルの「P＝」部分を「P＜」に改変する。



第4部　種間変異の解析

生物の種間変異とは、異なる種や亜種間で生じる遺伝的な差異のことで、生物の生息環境の違いもその差に寄与することがある。二枚貝はその生活様式を大まかに次の3種類に分類することが出来る：①堆積物に潜るタイプ、②堆積物の表面で生活するタイプ、③基盤に付着するタイプ。ここでは、①の生活様式を持つ種をAグループ、②と③の生活様式を持つ種をBグループとして、AとB二つのグループ間に分類し、形態の種間変異についての解析を行う事で、生活様式と形態の関係性について考察を行う。

解析用まとめデータセットの作成 動画8

まず、自分のデータの、「変量の比と代表値」のシートから、L/H、T/Hをコピーし、「解析用まとめデータ」シートの対応する場所にペーストする。



次に、「変量の比と代表値」のシートで左右の殻厚の違いの指標として、(TR-TL)÷(TRとTLの平均値)を計算して、コピーする。



コピーしたデータを、「解析用まとめデータ」シートの「(TR-TL)/mean」に値で貼り付ける。



「解析用まとめデータ」シートのHabitat_typeに適切な生活様式（AまたはB）を記入する。Sp_codeは種名コードの意味で、属名と種小名の頭文字を合わせて作成し、記入する。例えば、Patinopecten yessoensisであれば、Pyとする。



次に、WebClassの「解析用データ」から、リンクをたどり、以下の6種分のデータをダウンロードする。

＊解析用データ A_Anadara granosa.xlsx、A_Meretrix meretrix.xlsx、A_Scapharca subcrerate.xlsx、B_Argopecten irradians.xlsx、B_Chlamys larvata.xlsx、B_Comptopallium radula.xlsx

A_…で始まるデータは生活様式Aグループ、B_…で始まるデータは生活様式Bグループに属する種のものである。自分が計測した種がAグループに属する場合は、自分で測定したデータに加えて、自分が測定したものとは異なる2種分のAグループのデータとBグループ3種分のデータを用いる。自分が計測した種がBグループに属する場合は、自分で測定したデータに加えて、自分が測定したものとは異なる2種分のBグループのデータとAグループ3種分のデータを用いる。いずれの場合でも、合計6種分のデータを用いて以下の解析を行う。

各種のデータを、「解析用まとめデータ」シートの自分のデータの下にコピー＆ペーストする。全てのデータを貼り付けたら、「(TR-TL)/mean」となっている列名をDA（＊Directional Asymmetry、方向性のある左右非対称性の意味）に置き換えて、データセットを完成させる。



4-1. 分散分析と多重比較法

4-1-1. 分散分析

分散分析（ANOVA） は、複数のグループの 平均値に差があるか を調べる統計手法で、特に、3つ以上のグループを比較する場合によく使われる。

この方法では、データをグループごとに分類し、以下の2つのばらつきを比較する。

  • グループ内の変動：同じグループ内の個体のばらつき
  • グループ間の変動：グループ間の平均値の差

ばらつきの差が大きいほど、平均値の差が統計的に有意である可能性が高くなる。

手順は以下のとおり：

  1. 全変動（全体のばらつき）を求める
  2. それを グループ内変動 と グループ間変動 に分ける
  3. 両者の比から F値 を算出する
  4. F値を用いて、帰無仮説（グループ間の平均値に差がない）を検定する

4-1-2. 多重比較とTukey検定

多重比較法 は、複数のグループをペアで比較する際に用いられる。このとき問題になるのが 第1種の過誤（Type I error）、すなわち「誤って帰無仮説を棄却してしまうこと」。

そのリスクを抑える方法のひとつが Tukey検定。これは、分散分析の後に行い、各グループ間の平均値の差が有意かどうかを調べるのに使われる。

4-1-3. アルファベットによる有意差の表現

Tukey検定などの結果は、次のように アルファベット記号 を使って表現することがある：

  • 同じアルファベット → 統計的に有意な差がない
  • 異なるアルファベット → 有意な差がある

この記法は、視覚的に結果をわかりやすく示すためのもの。

L/HとDAについて、6種間で平均値に違いがあるかどうかを、分散分析とTukey検定を用いて解析せよ。また、平均値±標準誤差の棒グラフを作成し、種間のペアワイズ比較の結果をアルファベットで示せ（課題10、Excelファイル）。さらに、その結果の生態学的な解釈を記述せよ（課題11、Wordファイル）。

分散分析 動画9

これまでと同様に、「解析用まとめデータ」シートのデータ行列を格納する。データの読み込みには、以下のコードを用いる。

data <- read.table("clipboard", header = TRUE)

＊Macの場合にはこちら：　data<-read.table(pipe(“pbpaste”),header=T)



種名と生活様式グループの対応が分かりやすくなるように、新しい変数を作成する。

group<-as.factor(paste(data$Habitat_type,"_",data$Sp_code))
data<-data.frame(data,group)

このデータセットを用いて分散分析を実行する。

#L/H(結果を例示)
model<-aov(L.H~group,data=data)
summary(model)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## group         5 2.8747  0.5749   335.9 <2e-16 ***
## Residuals   144 0.2465  0.0017                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#DA
model2<-aov(DA~group,data=data)
summary(model2)

得られた結果のp-valueの値をコピーして、「分散分析の結果」シートの該当箇所（「分散分析　P=」）の部分に貼り付ける。

棒グラフ（エラーバーつき）の作図（Rを用いるバージョン） 動画10

tapply関数を用いて、平均値、標準偏差、標本数を算出し、標準誤差（SE: Standard Error）を計算する。

#L/H
m_L.H<-tapply(data$L.H, list(data$group), mean) #平均値
sd_L.H<-tapply(data$L.H, list(data$group), sd) #標準偏差
N_L.H<-tapply(data$L.H, list(data$group), length) #標本数
SE_L.H<-sd_L.H/sqrt(N_L.H) #標準誤差
#DA
m_DA<-tapply(data$DA, list(data$group), mean) #平均値
sd_DA<-tapply(data$DA, list(data$group), sd) #標準偏差
N_DA<-tapply(data$DA, list(data$group), length) #標本数
SE_DA<-sd_DA/sqrt(N_DA) #標準誤差

次に、平均値とSEを用いてエラーバー付き棒グラフを作図する。

#L/H(結果を例示)
BP_L.H<-barplot(m_L.H, ylim=c(0, max(m_L.H)*1.5))
arrows(BP_L.H,m_L.H+SE_L.H,BP_L.H,m_L.H-SE_L.H,code=3,angle=90,lwd=1,length=0.1)
box(lty=1)

#DA
BP_DA<-barplot(m_DA, ylim=c(min(m_DA)*1.5, max(m_DA)*1.5))
arrows(BP_DA,m_DA+SE_DA,BP_DA,m_DA-SE_DA,code=3,angle=90,lwd=1,length=0.1)
box(lty=1)

コードを実行して得られた棒グラフを「分散分析の結果」シートの該当箇所に貼り付ける。

Tukey検定 動画11

Tukey検定を行うために、multcompパッケージのインストールを行う。

Rのウィンドウ上部の「パッケージ」から「パッケージのインストール」と進む。



「Secure CRAN mirrors」からJapan（Yonezawa）を選択し、「Packages」から「multcomp」を選して、インストールする。



このmulticompを用いて、Tukey検定を実行する。

library(multcomp)#使用するパッケージの読み込み
#L/H(結果を例示)
L.H_res<-glht(model,linfct=mcp(group="Tukey")) 
summary(L.H_res)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = L.H ~ group, data = data)
## 
## Linear Hypotheses:
##                       Estimate Std. Error t value Pr(>|t|)    
## A _ Py - A _ Mm == 0 -0.164416   0.011701 -14.051  < 1e-04 ***
## A _ Ss - A _ Mm == 0  0.129640   0.011701  11.079  < 1e-04 ***
## B _ Ai - A _ Mm == 0 -0.109892   0.011701  -9.391  < 1e-04 ***
## B _ Cl - A _ Mm == 0 -0.259144   0.011701 -22.147  < 1e-04 ***
## B _ Cr - A _ Mm == 0 -0.252336   0.011701 -21.565  < 1e-04 ***
## A _ Ss - A _ Py == 0  0.294056   0.011701  25.130  < 1e-04 ***
## B _ Ai - A _ Py == 0  0.054524   0.011701   4.660 0.000101 ***
## B _ Cl - A _ Py == 0 -0.094728   0.011701  -8.096  < 1e-04 ***
## B _ Cr - A _ Py == 0 -0.087920   0.011701  -7.514  < 1e-04 ***
## B _ Ai - A _ Ss == 0 -0.239532   0.011701 -20.471  < 1e-04 ***
## B _ Cl - A _ Ss == 0 -0.388784   0.011701 -33.226  < 1e-04 ***
## B _ Cr - A _ Ss == 0 -0.381976   0.011701 -32.644  < 1e-04 ***
## B _ Cl - B _ Ai == 0 -0.149252   0.011701 -12.755  < 1e-04 ***
## B _ Cr - B _ Ai == 0 -0.142444   0.011701 -12.173  < 1e-04 ***
## B _ Cr - B _ Cl == 0  0.006808   0.011701   0.582 0.992094    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
cld(L.H_res)
## A _ Mm A _ Py A _ Ss B _ Ai B _ Cl B _ Cr 
##    "a"    "b"    "c"    "d"    "e"    "e"
#DA
DA_res<-glht(model2,linfct=mcp(group="Tukey")) 
summary(DA_res)
cld(DA_res)

ペアワイズ比較の結果のアルファベット表記が得られる。

これを「分散分析の結果」シートの対応する棒グラフにテキストファイルで書き込む。



[追加課題3、Excelファイル] 上記の解析コードを適宜改変して、T/Hについても同様の解析・作図を行ってみること。

＊生物科学実験IIの履修者は必須の課題とする。

[追加課題4、Excelファイル] 本実習のデータを用いて、本実習で描いたグラフ（ヒストグラム、棒グラフ）以外のグラフを２つ自由に描くこと。Excelを用いてもRを用いても良いが、描画結果はExcelファイルに貼り付けること。また、縦横軸のラベルを忘れずにつけること。

＊生物科学実験IIの履修者は必須の課題とする。